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Abstract — This  paper  presents  a  statistical-mechanics-inspired 
procedure  for  optimization  of  the  sensor  field  configuration  to 
detect  mobile  targets.  The  key  idea  is  to  capture  the  low¬ 
dimensional  behavior  of  the  sensor  field  configurations  across  the 
Pareto  front  in  a  multiobjective  scenario  for  optimal  sensor  de¬ 
ployment,  where  the  nondominated  points  are  concentrated  within 
a  small  region  of  the  large-dimensional  decision  space.  The  sensor 
distribution  is  constructed  using  location-dependent  energy-like 
functions  and  intensive  temperature-like  parameters  in  the  sense 
of  statistical  mechanics.  This  low-dimensional  representation  is 
shown  to  permit  rapid  optimization  of  the  sensor  field  distribu¬ 
tion  on  a  high-fidelity  simulation  test  bed  of  distributed  sensor 
networks. 

Index  Terms — Gibbs  distribution,  mobile  target  detection,  opti¬ 
mization  of  the  sensor  field  configuration,  sensor  networks. 

Nomenclature 

ek(i)  Scalar  energy  corresponding  to  the  kth  intensive  pa¬ 
rameter  and  the  ith  state. 

Ek  Energy  vector  corresponding  to  the  kth  intensive 
parameter. 

fs  Sensor  density  function. 

fr  Distribution  of  expected  targets. 

i  Index  for  the  state  (i  E  {1,2,...,  n}). 

J  Combined  objective  functional. 

k  Index  for  the  intensive  parameters  (k  E  { 1 ,  2 , . . . ,  K }). 

K  Total  number  of  intensive  parameters. 

kd  Number  of  sensors  confirming  target  detection. 

n  Number  of  states. 

N  Number  of  sensors. 
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Pd  Probability  of  detection  of  a  target  by  a  single  sensor 
(target  within  distance  Rd  from  the  sensor). 

Pp  Probability  of  a  single  sensor  false  alarm  rate. 

Pft  Probability  of  a  specific  false  track. 

Pfs  Probability  of  a  false  search. 

Pst  Probability  of  detection  of  a  specific  track. 

Pss  Probability  of  a  successful  search. 

Rd  Radius  of  detection  of  a  sensor. 

S  Surveillance  region. 

Wi  Gaussian-mixture  weight  corresponding  to  state  i. 

W  Weights  in  vector  form. 

W*(a)  Optimal  weights/distribution  for  a  given  a. 

Xi  x-position  of  the  Gaussian  mean  for  state  i. 

yi  ^/-position  of  the  Gaussian  mean  for  state  i. 

a  Tradeoff  variable  between  log  (Pss)  and  log  ( Pp  s )  • 

f3k  kth  intensive  parameter. 

( t>nT  Probability  of  finding  sensors  within  the  region 

VLT(x,y,6). 

A i  ith  largest  eigenvalue  of  H. 

rj  Threshold  for  approximation. 

Qt  Pill-shaped  region  around  a  target. 
a  Variance  of  each  component  in  the  Gaussian  mixture. 

I.  Introduction 

RECENT  advances  in  sensor  technology  and  the  develop¬ 
ment  of  powerful  mobile  computational  platforms  have  led 
to  the  usage  of  distributed  sensor  networks  for  the  tracking  of 
moving  targets  (e.g.,  undersea  autonomous  vehicles  and  wea¬ 
pon  systems)  [1]— [3].  Typically,  such  sensor  networks  consist  of 
a  large  number  (e.g.,  ^100-1000)  of  inexpensive  sensor  nodes 
that  cover  large  surveillance  regions.  As  an  advantage  of  the  par¬ 
allel  operations,  distributed  sensor  networks  have  a  larger  cov¬ 
erage  than  their  conventional  counterparts  in  terms  of  both  area 
and  speed.  In  addition,  distributed  sensor  networks  are  robust 
to  random  variations  in  operating  conditions  and  unexpected 
failures  that  may  happen  during  the  course  of  operation  [2]. 

Typically,  in  a  sensor  field,  each  sensor  is  designed  with 
a  limited  autonomous  target  detection  capability  to  reduce 
the  sensor  cost,  and  the  information  from  multiple  sensors  is 
combined  to  compensate  for  the  limited  coverage  provided 
by  the  individual  sensors.  This  process  of  requiring  multiple 
sensor  signals  for  a  successful  target  detection  is  referred  to  as 
track-before-detect  [4],  since  a  confident  system-level  detection 
decision  is  not  made  until  a  set  of  multiple  sensor  detections 
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occur  in  a  spatiotemporal  pattern  that  is  consistent  with  the  ex¬ 
pected  target  motion.  Along  this  line,  Wettergren  [5]  analyzed 
the  performance  of  the  track-before-detect  schemes  for  the 
sensor  networks.  Furthermore,  a  method  that  is  used  to  optimize 
the  sensor  field  configuration  for  an  efficient  target  detection 
was  proposed  using  different  performance  metrics  [6]. 

A  particular  problem  of  interest  is  that  of  the  optimization 
of  the  sensor  field  configuration  for  the  purpose  of  target  track 
coverage.  In  contrast  to  point  coverage  that  requires  optimal 
search  solutions  in  finding  stationary  objects  within  the  sur¬ 
veillance  region,  target  track  coverage  uses  fixed  sensors  to 
search  for  moving  targets  in  the  surveillance  region.  The  sensor 
density  that  is  required  for  target  track  coverage  is  lower  than 
that  for  point  coverage  as  the  moving  targets  leave  a  larger 
footprint  over  time,  and  therefore,  they  may  be  detected  by 
multiple  sensors.  Stone  [7]  presented  a  survey  of  optimal  target 
search  techniques,  where  only  the  independent  noncollabora- 
tive  sensors  were  used.  Meguerdichian  et  al.  [8]  examined  the 
performance  of  an  ad-hoc  network  to  address  the  problem  of 
sensor  allocation  to  achieve  the  target  tracking  capabilities  over 
the  region  of  interest.  Ram  et  al.  [9]  analyzed  the  track-coverage 
property  of  a  random  sensor  network  and  provided  asymptotic 
results  for  its  performance.  Cloqueur  et  al.  [10]  incorporated  a 
deployment  cost  to  choose  the  number  of  sensors  in  a  region. 
Track-coverage  assessment  of  a  given  configuration  was  ad¬ 
dressed  by  Baumgartner  and  Ferrari  [1 1]  for  the  reorganization 
of  the  sensor  field  to  achieve  the  maximum  coverage.  The  track- 
before-detect-based  optimal  control  of  sensor  reorganization 
was  analyzed  by  Baumgartner  et  al.  [12].  Recently,  the  concepts 
of  statistical  mechanics  [13]  have  been  extended  to  various 
fields  such  as  statistical  learning  [14],  communication  networks 
[15],  autonomous  systems  [16],  [17],  swarm  control  [18],  [19], 
and  optimization  [20],  [21]. 

Optimal  configurations  work  well  for  problems  of  limited 
duration.  However,  for  long-term  operation,  the  tradeoff  be¬ 
tween  detection  and  false  alarm  performance  may  change 
throughout  the  deployment.  In  such  cases,  the  reorganization 
of  these  fields  is  important  for  persistent  operation  to  accom¬ 
modate  changes  in  the  operational  intent.  The  reorganization 
of  the  sensor  field  requires  reoptimization  to  determine  its  new 
configuration.  However,  in  computation  and  communication 
constrained  environments,  computationally  efficient  methods 
for  rapid  optimization  are  required.  In  this  regard,  this  paper 
presents  a  statistical-mechanics-inspired  method  for  the  rapid 
optimization  of  the  sensor  field  configuration  that  provides  a 
set  point  for  the  new  configuration.  A  major  contribution  of 
this  paper  is  the  development  of  an  analytical  tool  to  facilitate 
the  rapid  multiobjective  optimization  of  the  sensor  placement 
for  the  detection  of  mobile  targets.  In  this  regard,  the  large¬ 
dimensional  space  of  the  sensor  field  configurations  is  ex¬ 
pressed  as  a  function  of  a  low-dimensional  set  of  intensive 
parameters  (e.g.,  temperature,  pressure,  and  chemical  poten¬ 
tial  in  statistical  mechanics).  Consequently,  the  optimization 
process  becomes  rapid  because  the  search  for  the  optimal  con¬ 
figuration  is  performed  in  a  significantly  reduced  dimensional 
space.  The  algorithms  for  the  optimization  of  the  sensor  field 
configuration  have  been  validated  on  a  simulation  test  bed  of 
distributed  sensor  networks. 


Fig.  1.  Detection  region  Qt(x,i/,6)  around  a  target  track  originating  at 
(x,y)  and  heading  along  a  direction  0  for  a  time  interval  St. 

This  paper  is  organized  in  six  sections  (including  the  current 
section)  and  an  Appendix.  Section  II  describes  the  detection 
model  that  is  used  to  obtain  the  probabilities  of  the  success¬ 
ful  and  false  search.  Section  III  poses  the  problem  of  rapid 
optimization  based  on  the  concepts  of  statistical  mechanics. 
Section  IV  formulates  the  optimal  algorithms  for  the  sensor 
field  configuration  in  the  setting  of  the  Gibbs  distribution,  and 
it  is  supported  by  the  Appendix.  Pertinent  results  are  presented 
in  Section  V.  This  paper  is  concluded  in  Section  VI,  with 
recommendations  for  future  research. 

II.  Performance  Model  for  Target  Detection 

This  section  describes  a  model  [5]  for  the  target  detection 
performance  of  a  typical  distributed  sensor  network.  The  global 
performance  is  characterized  by  the  following  two  measures: 

1)  the  probability  of  a  successful  search  (Pss),  which  is 
defined  by  the  multisensor  target  detection; 

2)  the  probability  of  a  false  search  (Pfs),  which  is  defined 
by  the  multisensor  false  alarms. 

While  the  details  of  computing  the  aforementioned  parame¬ 
ters  in  a  sensor  network  are  reported  in  the  previous  publi¬ 
cation  [5],  the  underlying  model  is  briefly  described  here  for 
completeness. 

A  collection  of  N  homogenous  sensors  is  deployed  in  a  given 
surveillance  region  <S,  where  it  may  be  necessary  to  update 
the  original  configuration  in  real  time  for  the  improvement 
of  the  detection  performance.  The  sensor  field  is  assumed  to 
have  the  following  characteristics:  1)  each  sensor  has  a  sensing 
radius  Ry,  within  which  it  may  successfully  detect  a  mobile 
target  with  a  probability  Pd  ,  and  2)  each  sensor  has  a  specified 
false  alarm  rate  (Pp),  i.e.,  the  probability  of  false  alarm  per 
unit  time.  The  probability  of  a  false  search  (Pfs)  is  reduced 
by  the  requirement  of  multiple  detections  by  different  sensors, 
occurring  in  a  sequence  that  is  spatiotemporally  consistent 
with  the  expected  target  motion  before  confirming  a  target.  In 
accordance  with  the  track-before-detect  paradigm  [4],  a  moving 
target  is  detected  if  the  ky  (typically  ky  =  3  or  4)  sensors  detect 
the  target  within  a  specified  time  interval  St  and  a  spatial  region, 
as  shown  in  Fig.  1.  This  time  interval  St  is  chosen  such  that 
the  expected  target  velocity  is  approximately  constant  over  the 
interval.  Such  an  approximation  is  consistent  with  the  Markov 
modeling  assumptions  of  the  moving  targets  that  are  commonly 
used  in  the  target  tracking  community  [22]. 
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Fig.  1  shows  a  target  track  originating  at  location  (x,  y )  and 
moving  at  a  heading  0  over  a  time  interval  St.  The  region 
Qt(x:  V ,  0)  around  this  target  track  corresponds  to  the  region  of 
target  detectability  over  time  interval  St.  This  region  represents 
a  subset  of  the  space  in  which  the  sensors  have  an  opportunity 
to  detect  the  specified  target  track  within  the  sensing  radius  Rd. 

Let  the  sensor  distribution  be  denoted  with  fs(x,y).  Then, 
the  probability  of  finding  a  sensor  in  a  region  Qt  is  given  by 

</>nT  =  J  fs(x',y')dx'dy'.  (1) 

(x',y')eQT 

Furthermore,  the  probability  of  a  single  sensor  detecting  the 
target  is  given  by  Pd&Qt  •  For  a  successful  target  detection,  it  is 
required  that  at  least  kd  out  of  N  sensors  independently  detects 
the  target  within  the  region  (each  with  equal  probability 
Pd^t\  This  is  modeled  by  the  Poisson  distribution  [2],  [5] 
for  a  large  numbers  of  sensors  N  and  is  given  as 

pST(ki,  a,) .  i  -  eM-NPoM  £ 

f  ml 

m= 0 

(2) 


Surveillance  region  S 


Intensive  parameters 
P1,  p2,  P3,  ...,PK 
common  to  all  states 


A 


State  i 
Weight:  W, 
Energy  parameters: 
e1  (i),  e2(i),...,  eK(i) 


(*i!  Vi) 

(X2.V2) 

(x3,'v3) 

(xi,  Vi) 

i 

y n) 

Fig.  2.  Partition  of  the  surveillance  region  and  the  associated  parameters. 

The  exact  procedure  for  the  evaluation  of  the  performance 
measures,  i.e.,  the  probability  of  a  successful  search  (Pss)  and 
the  probability  of  a  false  search  (Pps),  is  given  in  [5].  The 
procedure  for  determining  the  optimum  configuration  of  the 
sensor  field  for  detecting  the  moving  targets  is  described  in  [6]. 


III.  Optimization  Problem  for  Sensor  Placement 


where  PsT(kd^r)  is  the  probability  of  detecting  a  single 
target  that  moved  along  the  target  path  in  the  region  Q,T. 

Let  the  target  track  distribution  be  denoted  by  a  joint  prob¬ 
ability  density  function  fT(x,y,6),  i.e.,  the  likelihood  of  a 
target  being  present  at  location  (x,y),  with  heading  6.  Thus,  the 
overall  probability  of  a  successful  search  (Pss)  that  is  relative 
to  this  distribution  of  targets  is  given  by 

27 T 

Pss(kd)  =  f  f  PsT{kd,£lT)fTdxdydO.  (3) 

0  x,yeS 

Along  similar  lines,  the  probability  of  a  false  search  ( Pps ) 
is  defined  as  the  probability  that  at  least  kd  sensor  false  alarms 
occur,  which  are  kinematically  consistent  with  a  potential  target 
motion.  The  probability  of  a  single  sensor  raising  a  false  alarm 
in  a  time  interval  St  is  PpSt  for  the  small  values  of  the  false 
alarm  rate  Pp.  Similar  to  (2),  the  probability  of  a  specific  false 
track  in  the  region  Qt,  denoted  as  £2t),  is  given  by 


PFT(kd,  Ur)  =  1  -  exp(-NPF5t(l)nT ) 


(NPFSt^T 

2^  m! 


(4) 


The  probabilities  of  false  track  events,  which  are  independent 
for  the  nonoverlapping  regions  £2T,  are  assumed  to  be  Poisson. 
Thus,  the  probability  of  zero  false  tracks  is  given  by 


Pr{0}  =  exp 


■£//  PFT(kd,Q,T)dxdydo\  (5) 


0  x,yeS 


where  Aqt  is  the  area  of  the  region  QT.  Thus,  the  probability 
of  a  false  search  Pps  is  given  by 


PFs(kd)  =  l-Pr{0}.  (6) 


As  described  in  the  previous  section,  the  performance  of 
a  sensor  field  is  measured  in  terms  of  the  probability  of  the 
detection  of  the  moving  targets  and  the  associated  probability  of 
false  alarms.  The  false  alarm  objective  is  included  to  facilitate 
the  avoidance  of  sensor  clumping ,  which  may  cause  confusion 
in  a  track-before-detect  construct.  Additional  parameters  such 
as  tracking  accuracy  and  detection  time  may  also  be  used  in  the 
objective  functional.  The  optimization  of  the  sensor  distribution 
using  multiple  conflicting  objectives  leads  to  a  Pareto  optimal 
set ,  also  known  as  the  nondominated  set  [23]  of  the  sensor  field 
configurations.1 

Based  on  the  aforementioned  performance  measures  [see  (3) 
and  (6)],  the  objective  is  to  construct  a  set  of  nondominated 
sensor  distributions  fs(x,  y )•  In  this  regard,  these  distributions 
are  modeled  as  a  mixture  of  multivariate  Gaussian  density 
functions,  given  as 

{(x  -  Xi)2  +  (y-  Vi)2 
2  a2 

(7) 

where  n  is  the  number  of  components  in  the  Gaussian  mixture 
and  kFi,  VL2, . . . ,  Wn  are  the  mixture  weights.  The  Gaussian 
mixture  centers  (xi,yi)  are  chosen  a  priori  to  lie  on  a  square 
grid  on  the  search  space  S  consisting  of  n  cells,  as  shown 
in  Fig.  2.  Each  grid  cell  “i”  is  assigned  a  weight  W{.  The 
parameter  a  is  chosen  by  taking  into  account  the  grid  spac¬ 
ing  [6].  The  weights  are  normalized  to  yield  a  probability 
mass  function,  i.e.,  Y^i=i  Wi  =  1.  The  sensor  distribution 
fs  (x,  y)  can  be  obtained  by  computing  the  weight  vector  W  = 
[Wi,  W2, . . . ,  Wn\.  Therefore,  a  point  in  the  Pareto  set  repre¬ 
sents  a  particular  choice  of  the  weight  vector  W  as  used  in  (7). 

'The  configurations  in  the  Pareto  optimal  set  are  nondominated  in  the  sense 
that  there  is  no  configuration  which  is  an  improvement  in  every  objective  over 
any  Pareto  configuration. 


fs(x,y ) = 


27 RT2 


y:  Wi  exp 


2=1 
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Usually,  the  Pareto  set  is  a  very  small  subset  of  achievable 
configurations.  In  other  words,  a  vast  majority  of  the  sensor 
field  configurations  belong  to  the  dominated  set.  The  total  num¬ 
ber  of  possible  configurations  exponentially  increases  with  the 
number  of  sensors,  implying  that  the  Pareto  optimal  set  forms  a 
small  fraction  of  the  entire  set  of  configurations.  As  a  result, 
most  of  the  current  Pareto  optimization  techniques  [24]  that 
utilize  the  concepts  of  genetic  algorithms  or  nonlinear  program¬ 
ming  [23]  to  optimize  the  sensor  field  configuration  may  take 
a  long  time  to  search  the  entire  configuration  space  and  may 
still  yield  an  approximate  Pareto  set.  For  example,  traditional 
normal  boundary  intersection  (NBI)-based  genetic  algorithms, 
such  as  the  genetic-algorithm-based  NBI  (GANBI),2  [24]  may 
be  used  to  obtain  a  set  of  sensor  distributions  that  represent  a 
suboptimal  Pareto  set. 

This  paper  utilizes  the  solutions  obtained  from  GANBI  to 
create  a  model  of  the  Pareto  set  that  is  dependent  on  a  small 
number  of  parameters.  Consequently,  the  reduced  dimension 
of  the  search  space  enables  rapid  optimization  using  nonlinear 
programming  algorithms  (e.g.,  sequential  quadratic  program¬ 
ming  (SQP)  [23])  to  obtain  a  close-to-optimal  Pareto  set  which 
is  likely  to  be  superior  to  that  obtained  from  GANBI. 

To  define  the  objective  space,  this  paper  makes  use  of 
two  competing  performance  measures:  1)  the  probability  of  a 
successful  search  Pss{ W)  and  2)  the  probability  of  a  false 
search  Pfs( W),  where  W  denotes  the  weight  vector  that 
parameterizes  the  sensor  distribution  fs(x,y)  [see  (7)].  The 
other  parameters,  such  as  kd,  N,  and  fr{x,  V ,  0),  are  assumed 
to  be  implicitly  present  in  the  performance  measures.  Assuming 
convexity,  the  Pareto  front  is  parameterized  by  constructing  a 
maximization  objective  function  J  that  is  given  by  the  follow¬ 
ing  weighted  sum: 

J(W,  a)  =  a  log  [Pss( W)]  +(!-«)(-  log  [PFS( W)]) 


(8) 

where  aG  [0, 1],  i.e.,  the  Pareto  front  is  generated  by  varying 
a  from  zero  to  one.  The  Pareto  front  is  a  (generally  convex) 
curve  in  the  log (Pss)  versus  log (Pfs)  space.  In  this  space,  for 
a  point  on  the  Pareto  front  with  slope  m,  the  corresponding  a  = 
(1/1  -fra).  An  optimal  sensor  field  configuration  is  chosen 
from  the  Pareto  set  by  assigning  different  weights  to  different 
performance  objectives  [i.e.,  a  in  (8)].  Formally,  the  optimal 
sensor  distribution  for  a  given  value  of  a  is  given  by 

W*  {a)  =  arg  max  J( W,  a) .  (9) 

w 

In  real-time  operations,  a  supervisor  may  change  the  relative 
weights  of  the  individual  objectives  depending  on  the  situation, 
thereby  selecting  a  different  point  on  the  Pareto  set.  Conse¬ 
quently,  a  different  sensor  field  configuration  is  generated. 

Remark  3.1:  The  objective  function  J  in  (8)  may  be  aug¬ 
mented  with  additional  performance  measures,  such  as  tracking 


2GANBI  produces  a  set  of  nondominated  solutions  in  the  multiobjective 
space  as  an  estimate  of  the  Pareto  front. 


accuracy  and  detection  time.  However,  with  more  than  two 
objective  functions,  a  will  become  a  vector  instead  of  being 
a  scalar.  Correspondingly,  the  dimension  of  the  vector  a  will  be 
one  less  than  the  number  of  individual  objectives,  resulting  in 
a  Pareto  hypersurface.  The  rest  of  the  analysis  described  in  this 
paper  remains  unchanged. 

This  paper  utilizes  the  concepts  of  statistical  mechanics  in 
constructing  a  reduced-order  model  of  the  Pareto  set  with 
a  small  number  of  parameters.  This  model  provides  a  more 
tractable  structure  for  rapid  optimization  using  nonlinear  pro¬ 
gramming  methods.  Therefore,  in  addition  to  generating  a 
close-to-optimal  Pareto  front,  this  approach  enables  rapid  es¬ 
timation  of  the  parameters  for  sensor  field  configuration  when¬ 
ever  the  relative  weight  a  of  the  objectives  is  altered. 

The  underlying  principle  of  statistical  mechanics  involves 
the  construction  of  the  energy  states,  where  the  equilibrium 
probability  distribution  is  estimated  by  maximizing  the  entropy 
of  the  system  for  a  given  macroscopic  parameter  (e.g.,  energy) 
[13].  The  distribution  of  the  population  among  various  energy 
states,  called  as  the  generalized  canonical  distribution  or  the 
Gibbs  distribution  [13],  is  given  as 


exp  [~(3e(i)\ 
£"=iexp[-/3e(j)]’ 


ie  {l,2,...,n}  (10) 


where  e(i)  is  an  extensive  parameter  that  represents  the  energy 
of  a  particle  in  state  i  and  (3  is  an  intensive  parameter  that 
represents  the  inverse  temperature.  In  general,  a  system  is 
characterized  by  multiple  intensive  parameters  [13].  In  that 
case,  (10)  is  modified  as 


Pi  = 


exp 


-ZLi  Pkek(i) 


£j=i  exP  ~Ek=iPkekti) 


G  {1,  2, . . . ,  n} 

(ID 


where  (3k ,  with  k  G  {1,  2, . . . ,  K},  is  the  intensive  parameters 
and  ek(i ),  with  k  G  {1,  2, . . . ,  K},  is  the  extensive  parameters 
that  represent  the  energies  of  the  state  i. 

In  this  paper,  the  Gibbs  distribution  is  used  to  construct  a 
model  structure  for  the  rapid  optimization  of  the  sensor  field 
configurations.  The  optimization  procedure  consists  of  two 
phases:  1)  (offline)  training  and  2)  (online)  operation.  For  a 
sensor  field,  the  search  space  S  is  partitioned  into  n  cells  to 
form  a  grid,  and  each  grid  cell  is  defined  as  a  state,  as  shown  in 
Fig.  2.  In  the  training  phase,  the  K  energy  parameters  ek(i ), 
with  k  G  {1,  2, . . . ,  AT},  are  computed  for  each  state  i  (see 
Section  IV  for  details).  Once  computed,  the  energy  parameters 
are  fixed  for  each  state  and  are  kept  invariant  in  the  operation 
phase.  The  intensive  parameters  (3k ,  with  k  G  {1,2 , . . . ,  if}, 
do  not  depend  on  the  state  i.  However,  they  may  vary  in  the 
operation  phase.  In  the  operation  phase,  moving  to  a  different 
operating  point  on  the  Pareto  front  corresponds  to  a  change  in 
the  intensive  parameters  while  keeping  the  associated  energies 
ek(i )  invariant.  Since  the  intensive  parameters  are  common  to 
all  states,  they  serve  as  the  control  parameters  for  the  optimiza¬ 
tion  of  the  sensor  field  configuration,  which  may  be  broadcasted 
to  the  entire  network  for  online  operation  [25]. 
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IV.  Gibbs  Distribution  for  Optimizing  the 
Sensor  Field  Configuration 

This  section  describes  a  method  that  is  used  to  model  the 
changes  in  the  sensor  distribution  as  the  parameter  a  in  (8) 
varies  in  the  interval  [0,  1].  It  is  desirable  to  construct  this 
model  of  the  Pareto  set  with  as  few  parameters  as  possible  in 
order  to  facilitate  the  real-time  adaptation  for  the  sensor  field 
configuration.  Any  point  on  the  convex  Pareto  front  can  be 
represented  as  an  optimal  solution  that  maximizes  J(W,  a) 
in  (8)  for  some  a  E  [0,1].  Therefore,  a  continuous  one-to-one 
mapping  exists  between  the  Pareto  front  and  the  set  [0,  1]. 
As  described  in  (7),  the  sensor  distribution  fs(x,y)  is  pa¬ 
rameterized  by  the  weights  W\ ,  W2 , . . . ,  Wn  given  to  the  n 
components  in  the  Gaussian  mixture.  Furthermore,  each  point 
on  the  Pareto  front  corresponds  to  a  certain  weight  vector 
expressed  as  W* (a)  =  [  {a)  W2*  (a)  •  •  •  W* {a) } . 

A.  Model  Construction 

The  individual  weights  are  modeled  by  the  Gibbs  distribu¬ 
tion  as 

W* (a)  =  exp  /3k(a)ek(i)  —  (12) 

where  the  energy  functions  ek(i),  with  k  E  {1, . . . ,  AT}  and 
K  <n,  are  defined  on  states  i  E  {1, . . . ,  n}  and  A(a)  is  the 
normalizing  factor.  Note  that  the  intensive  parameters  f3k ,  with 
k  E  {1, . . . ,  if},  depend  on  a  but  not  on  the  state  i,  while  the 
extensive  parameters  ek ,  with  k  E  {1, . . . ,  if },  depend  on  the 
state  i  but  not  on  a.  Taking  the  logarithm  on  both  sides  and 
multiplying  with  the  negative  one,  (12)  yields 

K 

Gi(a )  =  -  log  (W;(a))  =  £  P\a)ek(i)  +  A(a)  (13) 

k= 1 

where  (3k(a)  and  A  (a)  depend  on  the  parameter  a.  The  term 
A(a)  that  is  chosen  to  normalize  the  distribution  is  analogous 
to  the  free  energy  in  the  thermodynamic  sense  [13]. 

As  stated  earlier,  the  state-dependent  ek(i)’s  are  invariant 
along  the  Pareto  front,  while  (3k’s  are  invariant  across  all 
states.  The  next  task  is  to  estimate  the  values  of  the  state- 
dependent  ek,s  in  the  sense  of  minimum  variance.  Since  these 
functions  are  defined  over  a  discrete  and  finite  domain,  they 
can  be  represented  by  the  corresponding  energy  vectors:  Ek  = 
[  ek  (1)  ek(2)  •••  ek(n)  ]T  Vfc  E  {1, ...,  if }.  The  follow¬ 

ing  constraints  are  imposed  without  loss  of  generality. 

1)  The  Euclidean  norm  of  the  energy  vector  is  set  as 
\\Ek  1 1 2  =  1  Vfc  E  {1, . . . ,  if}.  The  scaling  factor  of  the 
energy  vector  is  absorbed  into  f3k . 

2)  Any  constant  added  to  the  energy  vectors  is  absorbed 
into  the  bias  function  A(a ),  and  it  does  not  affect  the 
distribution.  Hence,  the  sum  of  the  elements  of  the  energy 
vector  is  constrained  to  be  zero  by  subtracting  the  mean 
of  the  energy  from  each  element.  Thus 

n 

]Tefc(i)  =  o  Vke{i,...,K}. 

2=1 


3)  Given  K  linearly  independent  energy  vectors  E1  .E2. 
...,EK,  it  is  possible  to  construct  an  orthogonal  set 
of  if  vectors  that  span  the  same  space.  In  this  set  of 
basis  vectors,  the  coefficients  f3k  {a)  are  calculated  by  an 
appropriate  transformation.  Therefore,  it  is  assumed  that 
the  energy  vectors  are  mutually  orthogonal.  Thus 

(Ekl)T(Ek2)  =  5k1k2  h,k2e{l,...,K}  (15) 

where  Sk1k2  is  the  Kronecker  delta  function. 

Upon  imposing  the  aforementioned  constraints  and  summing 
(13)  over  all  n  states,  it  follows  that 

n  K  n  n 

£<3,  (a)  =  ^2,f3k(a)'Y^ek{i)  +  ^A(a)  (16) 

2=1  k=  1  2=1  2=1 

1  71 

=>  A(a)  =  -  Gi(a)  by  (14).  (17) 


By  setting  Gi(a)  =  Gi(a)  —  A(a ),  a  vector  is  defined  as 


G(<r)  —  Gi(o')  G2  (ci) 


Gn(a) 


1  T 


The  components  of  vector  G(a)  also  sum  up  to  zero,  i.e., 
the  vectors  G(a)  and  El,E2, . . . ,  EK  lie  in  the  same  (n  —  1)- 
dimensional  subspace  of  Mn. 


B.  Offline  Estimation  of  the  Energy  Vectors 

In  the  training  phase,  the  energy  vectors  are  estimated  from 
a  collection  of  points  in  the  Pareto  set  to  represent  the  non- 
dominated  sensor  field  distributions.  In  this  paper,  these  points 
are  generated  by  a  genetic-algorithm-based  multiobjective  op¬ 
timization  algorithm  (GANBI)  [24]. 

Let  the  training  data  set  consist  of  L  samples,  and  let  the 
n-dimensional  weight  vectors  that  correspond  to  these  sensor 
distributions  be  given  by  W1 ,  W2 , . . . ,  WL.  For  a  reliable  sta¬ 
tistical  analysis,  L  must  be  sufficiently  larger  than  the  number 
of  states  n.  Let  G^  be  a  transformation  of  the  weight  vector 
such  that 


Gl  =  —  log(Wz),  l  e  (18) 


As  defined  earlier,  the  vector  Gl  is  obtained  as 


W  A 


-  G'  —  \-EGi  )  [M,  •••)!] 


T 

nxl) 


2=1 


l  £  {1, ...  ,L}. 

(19) 


The  orthonormal  vectors  E1,  E2 , . . . ,  EK  are  estimated  by 
minimizing  the  sum  of  the  squared  distances  between  the 
vectors  Gl  and  their  projections  onto  the  space  spanned  by 
E1,  E2, . . . ,  Ek  .  To  this  end,  a  cost  functional  is  defined  as 


*  =  £ 


i=i 


K  _ 
G'  -  J2Ek  ((Ek)T G*) 


k=l 


(20) 


(14)  where  £  is  minimized  by  obtaining  its  vector  derivatives  with 
respect  to  the  energy  vectors  E1 ,  E2, . . . ,  EK  and  by  setting 
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them  to  be  equal  to  zero.  The  orthonormality  constraints  are  in¬ 
corporated  by  means  of  Lagrange  multipliers  (see  the  Appendix 
for  details).  The  energy  vectors  are  obtained  as  the  eigenvec¬ 
tors  of  the  self-adjoint  and  positive  semidefinite  matrix  H  = 

(Et  i(Gl)(Glf)  such  that 


HE  =  A  E. 


(21) 


The  vectors  E1 ,  E2 , . . . ,  EK  obtained  by  the  discrete 
Karhunen-Loeve  expansion  [26]  are  real  and  mutually  orthog¬ 
onal.  The  eigenvalues  are  ordered  asAi  >  A2  >•••  >  An  >0. 
The  eigenvectors  corresponding  to  the  largest  K  eigenvalues, 
where  K  <  n,  are  chosen  such  that  the  cumulative  sum  (77) 
of  the  remaining  (n  —  K)  eigenvalues  is  bounded  above  by  a 
positive  threshold  rjT  1 


A  Ysi=K+ 1  ^ 

v~  E"=1  Ai 


—  Vt- 


(22) 


As  described  earlier,  the  energy  vectors  corresponding  to 
each  state  form  the  local  potentials,  and  they  are  identified 
a  priori  in  the  (offline)  training  phase.  The  training  data  set  of 
the  nondominated  sensor  distributions  described  by  the  weight 
vectors  Wl  (l  =  1, . . . ,  L)  are  obtained  using  the  genetic  al¬ 
gorithms.  This  data  set  is  suboptimal  compared  to  the  ideal 
Pareto  front.  Although  the  values  of  f3k ,  with  k  =  1, . . . ,  K, 
may  be  found  for  each  sensor  distribution  in  the  training  set, 
the  weight  vectors  (l  =  1, . . . ,  L)  are  not  labeled  with  the 
corresponding  values  of  a,  making  it  impossible  to  obtain  (3k  as 
a  function  of  a.  This  issue  is  addressed  in  the  (online)  operation 
phase. 


C.  Online  Operation 


The  probability  of  a  sensor  occupying  a  state  i  is  governed 
by  a  combination  of  global  parameters  /3k,s  and  local  energy 
parameters  ek(i)’ s.  Given  the  energy  vectors,  as  obtained  in 
the  training  phase,  the  weight  vector  W  depends  on  the  values 
of  the  intensive  parameters  (3k ,  with  k  =  1, . . . ,  K,  and  it  is 
given  by 


exP  (-£fc=i  Pkek(i)) 

£”=iexp(-£f=1/3^(j))’ 


i  =  (23) 


where  K  is  much  smaller  than  the  number  of  states  n.  For 
a  given  value  of  a,  the  objective  functional  J(a)  [see  (8)]  is 
constructed  as  a  linear  combination  of  the  probabilities  of  the 
successful  and  false  searches.  The  aim  here  is  to  find  the  weight 
vector  W(<y)  that  maximizes  the  objective  J(a).  However,  the 
weight  vectors  are  functions  of  f3k ,  with  k  =  1, . . . ,  K,  that 
belong  to  an  abstract  AT-dimensional  space.  The  task  of  finding 
the  optimal  weight  vector  W  (a)  is  reduced  to  that  of  finding  the 
optimal  parameters  /3k(a),  which  maximize  J(a).  Due  to  the 
small  dimensionality  of  the  search  space  (i.e.,  the  space  of 
f3k’ s),  the  SQP  algorithm  is  suitable  in  computing  /3k(a).  Thus, 
the  supervisory  controller  may  choose  an  operating  point  in  the 
AT-dimensional  space  of  /3k’s  to  achieve  the  global  objectives 
for  a  given  a. 


V.  Results  and  Discussion 

This  section  presents  the  results  generated  by  using  the 
multiple-intensive  parameter  Gibbs  distribution  to  model  the 
Pareto  front  of  the  nondominated  sensor  field  configurations. 
The  competing  performance  criteria  in  the  Pareto  front  are  the 
probability  of  a  successful  search  (Pss)  and  the  probability  of  a 
false  search  (Pfs)  for  the  moving  targets.  While  the  details  are 
reported  in  [6],  a  brief  description  for  the  evaluation  of  Pss  and 
Pfs  for  a  given  sensor  field  configuration  has  been  described 
in  Section  II.  In  the  test  problem,  the  surveillance  region  is 
divided  into  a  20  x  20  square  lattice  to  yield  n  =  400  grid 
cells.  The  sensor  field  configuration,  expressed  as  a  probability 
density,  is  modeled  by  a  Gaussian  mixture  centered  at  each 
of  the  grid  points  [see  (7)].  The  distribution  of  the  sensors 
is  completely  parameterized  by  a  400-element  weight  vector 
(i.e.,  n  =  400),  resulting  in  399  independent  parameters  be¬ 
cause  of  the  constraint  XlSi  Wi  =  1.  The  global  performance 
parameters  Pss  and  Pfs  have  been  evaluated  based  on  the 
following:  1)  a  uniform  distribution  of  the  target  tracks,  i.e., 
/t(x,t/,  0)  =  (1/27 rL2),  and  2)  a  skewed  distribution  of  the 
target  tracks  /t(x,t/,  0)  =  (1  +  cos(6) /2ttL2)  in  the  domain 
of  the  L  x  L  surveillance  region.  In  these  simulation  scenarios, 
the  total  number  of  sensors  required  to  be  placed  is  chosen  to  be 
N  =  60,  with  the  cumulative  point  coverage  of  all  the  sensors 
amounting  to  only  9%  of  the  surveillance  region.  In  accordance 
with  the  track-before-detect  strategy,  the  number  of  sensors  kd 
that  must  detect  a  target  to  confirm  its  presence  is  chosen  as  two 
for  the  results  that  follow. 

A  GANBI  code  [24]  has  been  adopted  to  solve  the  multi¬ 
objective  optimization  problem  and  to  obtain  the  approximate 
Pareto  front  of  sensor  distributions,  where  8  b  are  used  to 
represent  each  weight  Wi  in  the  genetic  program  [6].  With 
a  population  size  of  100,  the  genetic  algorithm  has  been  run 
for  300  generations  to  obtain  a  good  spread  of  points  along 
the  Pareto  front,  where  each  point  corresponds  to  a  set  of 
400  weights  that  are  used  to  describe  the  sensor  distribution. 
As  with  any  genetic-algorithm-based  approach,  the  generated 
Pareto  front  is  always  approximate. 

The  algorithm  presented  in  this  paper  extracts  the  low¬ 
dimensional  characteristics  of  the  approximate  Pareto  front 
in  the  space  of  the  weights  Wi  and  formulates  a  tractable 
optimization  problem  with  fewer  parameters  in  order  to  refine 
and  improve  the  Pareto  front.  To  this  end,  the  energy  vectors 
E\E2,...,Ek  are  obtained  from  the  set  of  approximate 
Pareto  optimal  sensor  distributions  as  obtained  from  GANBI 
[24].  The  number  of  energy  vectors  K  is  given  by  (22)  after  an 
appropriate  choice  of  threshold  7]t- 

Once  the  energy  vectors  are  evaluated,  the  sensor  field  dis¬ 
tribution,  given  by  the  weights  Wi,  is  parameterized  in  terms 
of  the  intensive  parameters  /51 ,  /52 , . . . ,  /3K  [see  (23)].  Given 
an  objective  function  in  terms  of  a  [see  (8)]  and  keeping  the 
energy  vectors  invariant,  the  optimal  sensor  distribution  param¬ 
eters  Z?1,  ft2, . . . ,  /3K  are  obtained  by  using  the  SQP  algorithm 
as  follows: 

(f31(a)  ,...,pK(a))*  =  arg  max^  J  (W(/3\ . . . ,  pK),  a) . 

(24) 
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Fig.  3.  Sensor  density  per  unit  area  corresponding  to  the  points  on  the  Pareto 
front  (K  =  8)  for  the  uniform  target  distribution  fF(x,y,  6)  =  ( 1/2ttL 2). 
(a)  Configuration  A:  Pss  =  0.584  and  Pps  =  0.043.  (b)  Configuration  B: 
Pss  =  0.577  and  Pps  =  0.017.  (c)  Configuration  C:  Pss  =  0.549  and 
Pfs  =  0.011.  (d)  Configuration  D:  Pss  =  0.485  and  Pfs  =  0.007. 
(e)  Pareto  fronts  with  two,  six,  and  eight  intensive  parameters. 

A  complete  Pareto  front  is  generated  by  changing  the  scalar 
a  in  the  objective  function  [in  (8)]  over  a  desired  range  to 
obtain  the  optimal  values  of  the  parameters  /51 ,  /52 , . . . ,  /3K . 
The  bottom  plate  in  Fig.  3  shows  the  estimated  Pareto  fronts 
with  increasing  numbers  (If)  of  intensive  parameters  for  the 
uniform  target  distribution  fx(x,y,6)  =  (1/2ttL2).  Thus,  the 
399-dimensional  space  of  the  weight  vectors  is  reduced  to  a 
AT-dimensional  space  of  intensive  parameters.  The  sensor  den¬ 
sities  per  unit  area  displayed  in  the  four  plates  in  Fig.  3(a)-(d) 
correspond  to  the  four  points  marked  as  A,  B,  C,  and  D  on 
the  Pareto  front  in  Fig.  3(e),  respectively.  Similar  results  are 
shown  in  Fig.  4  for  a  skewed  target  distribution  fr(x,  y,  0)  = 
(1  +  cos(O) / 2tt L2) .  Although  the  Pareto  fronts  become  more 
accurate  with  a  larger  K ,  the  incremental  improvements  do  not 
justify  the  increased  computational  complexity  beyond  a  certain 
point  (in  this  case,  K  =  6).  The  errors  due  to  the  reduction  in 
the  order  of  the  decision  space  decrease  with  an  increase  in  the 
number  of  intensive  parameters,  as  shown  in  the  plots  in  Fig.  5, 
and  they  are  given  by  y  in  (22).  The  approximate  computation 


Fig.  4.  Sensor  density  per  unit  area  corresponding  to  the  points  on  the 
Pareto  front  (K  =  7)  for  the  skewed  target  distribution  fr(x,y,0)  = 
(1  +  cos(6>)/2ttL2).  (a)  Configuration  A:  Pss  =  0.581  and  PFS  =  0.035. 
(b)  Configuration  B:  Pss  =  0.575  and  PFs  =  0.023.  (c)  Configuration  C: 
Pss  =  0.544  and  PFs  =  0.016.  (d)  Configuration  D:  Pss  =  0.459  and 
PFS  =  0.010.  (e)  Pareto  fronts  with  three  and  seven  intensive  parameters. 


Fig.  5.  Error  due  to  a  reduced-order  decision  space. 

times  for  the  algorithms  are  shown  in  Table  I.  For  the  example 
described  in  Fig.  3,  the  genetic  algorithm  (GANBI)  used  to 
obtain  the  initial  Pareto  front  takes  ^6100  s,  and  it  is  usually 
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TABLE  I 

Comparison  of  the  Computation  Time  for  the  Algorithms 


GANBI 

Statistical  Mechanics-inspired 
Rapid  Optimization  Method 

K=2 

K=6 

K=8 

Optimization 
for  a  given  a 

N.Af 

3s 

20  s 

30  s 

Optimization 
for  Pareto  front 

6100  s 

60s6 

400  s6 

600  s6 

executed  offline.  Once  the  energy  vectors  are  obtained,  the 
procedure  for  optimizing  the  parameters  /31,  /32, . . . ,  /3K  is 
relatively  less  time  consuming.  Given  a  value  of  a ,  this  opti¬ 
mization  procedure  takes  2,  20,  and  30  s  for  K  —  2,  6,  and  8, 
respectively.  These  results  support  the  applicability  of  the  pro¬ 
posed  method  for  the  purpose  of  the  real-time  optimization  of 
the  sensor  placement  for  target  detection. 

Remark  5.1:  The  computation  time  for  the  execution  of  the 
GANBI  strongly  depends  on  the  number  of  states  n.  A  larger  n 
may  lead  to  a  greater  accuracy  in  determining  the  sensor  config¬ 
uration  at  the  cost  of  an  increased  computational  load  (i.e.,  exe¬ 
cution  time  and  memory  requirements).  On  the  other  hand,  the 
complexity  of  the  online  statistical-mechanics-inspired  method 
depends  on  the  number  of  intensive  parameters  K  that  can  be 
chosen  based  on  the  information  provided  in  Table  I. 

The  sensor  distributions  in  Figs.  3(a)  and  4(a)  correspond  to 
the  case  where  the  maximum  weight  is  given  to  the  probability 
of  a  successful  search  (i.e.,  a  =  1).  Fig.  3(a)  considers  the 
targets  that  may  move  in  any  direction  with  equal  probability. 
As  a  result,  the  sensor  distribution  is  circularly  symmetric 
and  is  in  the  shape  of  a  ring.  On  the  other  hand,  Fig.  4(a) 
considers  the  target  movements  predominantly  from  left  to 
right.  Consequently,  the  sensor  distribution  is  in  the  form  of 
barriers  spanning  from  top  to  bottom.  In  both  the  cases,  as 
the  weight  on  the  probability  of  a  successful  search  is  reduced 
(a  is  reduced),  the  sensor  distribution  distorts  are  a  manner  to 
effectively  reduce  the  probability  of  a  false  search. 

VI.  Summary  and  Conclusion 

This  paper  has  presented  a  methodology  for  the  rapid  op¬ 
timization  of  the  sensor  field  configuration  for  a  collaborative 
network  for  the  detection  of  the  moving  targets.  The  underlying 
algorithms  are  based  on  the  concepts  of  equilibrium  statistical 
mechanics,  where  the  Gibbs  distribution  has  been  used  to  model 
the  spatial  sensor  distribution.  In  particular,  the  number  of 
intensive  parameters  is  significantly  smaller  than  the  number  of 
states.  The  energy-like  quantities  are  spatially  dependent,  but 
they  are  independent  from  the  operating  point  on  the  Pareto 
front.  On  the  other  hand,  the  intensive  parameters  are  spatially 
invariant,  but  they  are  dependent  on  the  location  of  the  operat¬ 
ing  point  on  the  Pareto  front.  This  approach  to  suboptimal  rep¬ 
resentation  attempts  to  capture  the  low-dimensional  behavior 
of  the  nondominated  configurations  of  a  very  high  dimensional 
system.  The  following  topics  are  under  active  research. 

1)  The  sensitivity  analysis  of  the  sensor  distribution  with 
respect  to  the  intensive  parameters  and  its  effects  on  the 
computation  time  for  optimization. 


2)  The  extension  of  the  Gibbs  distribution  to  the  generalized 
exponential  distributions  for  modeling  the  Pareto  front. 

3)  The  analysis  of  the  computational  complexity  versus  the 
accuracy  tradeoff  in  the  optimization  process  by  choosing 
different  numbers  of  intensive  parameters. 

4)  The  reorganization  of  the  sensor  fields  in  real  time  to 
detect  the  mobile  targets  and  the  analysis  of  the  actual 
sensor  trajectories  during  the  transition.  The  reorgani¬ 
zation  may  also  be  addressed  by  sensor  activation  and 
deactivation  in  various  parts  of  the  field. 


Appendix 

Estimation  of  the  Energy  Vectors 


This  Appendix  outlines  a  procedure  for  estimating  the  energy 
vectors  Ek’s  by  minimizing  the  cost  function  £  defined  in  (20) 

l  /  k  \T 

£  =  E  G'-E^  ((£,fe)r G‘)  I 

1=1  \  k=l  ) 

x  {^l  ~EEk  {(Ek)TGl^j 

L  L  K 

=  J2(&)TGl  -2  ^2  '52(Gl)TEk(Ek)TGl 

1=1  1=1  k=l 

L  K  K 


+  EEE(G‘  )TEkl  ( Ekl  )TEk2  ( Ek 2  )TGl. 

1=1  fc1  =  1  k2  =  1 


The  usage  of  the  orthonormal  property  of  Ek ,  with  k  =  1 , . . . ,  K, 
yields 

£  =  f(G*)TG*  -  ^(Gz)T£fc(£*)TGz  j  . 

1=1  \  k=l  ) 

Since  the  quantities  (Gl)T  Ek  and  (Ek)TGl  are  scalars,  their 
order  of  multiplication  is  reversed  to  yield 

£  =  ^2  f(Gz)TG*  -  ^2(Ek)TGl(Gl)TE  )  • 

i=i  V  k= i  J 


A  Lagrange  function  £  is  constructed  by  introducing  new 
variables  A&,  i  =  1,  2, , . . ,  K,  to  incorporate  the  unit  length 
constraints  on  the  energy  vector  and  is  given  by 

K 

£  =  £  +  Y/^k((Ek)TEk  -!)■ 

k=l 


Taking  the  matrix  derivatives  of  £  with  respect  to  the  energy 
vectors  and  setting  them  equal  to  zero  yield 

^  =  -2 (Ekf  (^2(Gl)(Gl)T^  +  2A k(Ekf  =  0 

=>  Ek  =  A kEk ■ 

The  energy  vectors  are  estimated  as  the  eigenvectors  of  the 
matrix  (%2i= i{Gl)(Gl)T).  The  symmetry  property  of  the  ma¬ 
trix  implies  that  the  energy  vectors  are  orthogonal  to  each  other. 
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